import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用 SimHei 字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义微分方程组
def model(y, t, params):
    C, W, I, B, Bt = y  # 设置初始种群数量
    r_C, a_CI, f_C, g_C, d_C, r_W, d_W, e_IC, a_IC, m_IB, m_IBt, e_BI, a_IB, d_B, e_BtI, d_Bt = params

    # r_C:农作物的自然增长率            # a_CI:昆虫对农作物的攻击率
    # f_C:蝙蝠传粉使作物增长率           # g_C:杂草对作物的竞争率
    # d_C:农作物自然死亡率              # r_W:杂草的自然增长率
    # d_W:杂草自然死亡率                # e_IC:昆虫从农作物获得能量的效率
    # a_IC:昆虫捕食农作物的效率          # m_IB:鸟类捕食昆虫的效率
    # m_IBt:蝙蝠捕食昆虫的效率           # e_BI:鸟类从昆虫获得能量的效率
    # a_IB:鸟类捕食昆虫的效率            # d_B:鸟类的自然死亡率
    # e_BtI:蝙蝠从昆虫获得能量的效率     # d_Bt:蝙蝠的自然死亡率

    dCdt = r_C * C + f_C * Bt * C - a_CI * I * C - g_C * W * C - d_C * C  # 农作物增长
    dWdt = r_W * W - d_W * W - g_C * C * W  # 杂草增长
    dIdt = e_IC * a_IC * C * I - m_IB * B * I - m_IBt * Bt * I  # 昆虫增长
    dBdt = e_BI * a_IB * I * B - d_B * B  # 鸟类增长
    dBtdt = e_BtI * m_IBt * I * Bt - d_Bt * Bt  # 蝙蝠增长

    return [dCdt, dWdt, dIdt, dBdt, dBtdt]

# 初始种群数量
C0 = 9746  # 农作物
W0 = 324  # 杂草
I0 = 75  # 昆虫
B0 = 6  # 鸟类
Bt0 = 5  # 蝙蝠
y0 = [C0, W0, I0, B0, Bt0]

# 时间点
t = np.linspace(5, 15, 1000)

# 调整后的参数
params = (0.0005, 0.0003, 0.015, 0.0003, 0.065, 0.8, 0.01, 0.02, 0.0004, 0.01, 0.006, 0.1, 0.002, 0.005, 0.06, 0.015)

# 求解微分方程组
solution = odeint(model, y0, t, args=(params,))

# 绘制结果
plt.figure(figsize=(10, 8))
plt.plot(t, solution[:, 0] / C0, label='农作物数量 (C)')
plt.plot(t, solution[:, 1] / W0, label='杂草数量 (W)')
plt.plot(t, solution[:, 2] / I0, label='昆虫数量 (I)')
plt.plot(t, solution[:, 3] / B0, label='鸟类数量 (B)')
plt.plot(t, solution[:, 4] / Bt0, label='蝙蝠数量 (Bt)')

plt.legend()

plt.xlabel('时间')
plt.ylabel('当前种群数量/初始种群数量')
plt.title('农田生态系统成熟后种群规模随时间变化（停止使用化学物）')
plt.show()